import nglview
import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image,HTML
import __main__
from time import sleep
import pymol
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
На основе одного липида созадим ячейку с 64 липидами.
%%bash
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.reinitialize()
pymol.cmd.load('b_64.pdb')
pymol.cmd.do('png b_64.png')
Image(filename='b_64.png')
В файле b.top установите правильное количество липидов в системе (меняем 34 на 64)
! sed -i 's/34/64/g' b.top
Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды.
%%bash
editconf -f b_64.gro -o b_ec -d 0.5
Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул.
%%bash
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v
Запустив в консоли получил:
F{max}^0 = 4.37970e+05
F{max}^{1} = 6.16887e+02
Добавим в ячейку молекулы воды типа spc.
%%bash
genbox -cp b_em -p b -cs spc216 -o b_s
Проведём "утряску" воды:
%%bash
grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v
%%bash
grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.load('b_s.pdb')
pymol.cmd.do('png b_s.png')
Image(filename='b_s.png')
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.load('b_pr.pdb')
pymol.cmd.do('png b_pr.png')
Image(filename='b_pr.png')
Заметно, что вода более равномерно и плотно расположилась внутри ячейки. Сама ячейка стала "беспорядочна"
%%bash
grompp -f md -c b_pr -p b -o b_md -maxwarn 1
Анализ результатов вычислений на кластере¶
%%bash
echo 2 | trjconv -f b_md.gpu.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
import io
from base64 import b64encode
video = io.open('mo.mp4', 'r+b').read()
video_encoded = b64encode(video)
html_code = '<video controls alt="PyMol Movie" src="data:video/mp4;base64,{}" type="video/mp4">'.format(video_encoded)
display(HTML(html_code))
На 59 изображении отчетливо видно разделение на два слоя.
MODEL59 соответствует t=29000.00000
Image(filename='img59.png')
%%bash
echo 2 | g_traj -f b_md.gpu.xtc -s b_md.tpr -ob box_1.xvg
import pandas as pd
coord = pd.read_csv('box_1.xvg', header = None, skiprows = 24, sep='\t')
coord = coord[[1,2,3,4]]
coord.columns = ['t', 'x', 'y', 'z']
coord.head(10)
Площадь считаем как произведение по осям Y и Z и нормируем на 32 (число липидов в слое)
plt.figure(figsize=(15,8))
lip_sizes = np.array(coord.y)*np.array(coord.z)/32.
plt.plot(np.array(coord.t), lip_sizes);
Судя по графику со временем площадь липида уменьшается. В процессе моделирования структура компактно собирается, приходя в оптимальное состояние.
Посмотрим на площадь доступных гидрофильной и гидрофобной поверхностей.
%%bash
echo 2 | g_sas -f b_md.gpu.xtc -s b_md.tpr -o sas_b.xvg -dt 100
hydro = pd.read_csv('sas_b.xvg', header = None, skiprows = 22, sep='\s+')
hydro = hydro[[0,1,2]]
hydro.columns = ['t', 'Hydrophobic', 'Hydrophilic']
hydro.head(10)
plt.figure(figsize=(20,8))
plt.plot(np.array(hydro.t), hydro.Hydrophobic, label='Hydrophobic', color='red')
plt.plot(np.array(hydro.t), hydro.Hydrophilic, label='Hydrophilic', color='blue')
plt.legend()
Площадь доступной гидрофобной поверхности стремительно снижается в первые ~5000 шагов, что соответствует увиденному выше образованию гидрофобного ядра. Это логично: упаковка гидрофобной поверхности в ядро существенно снижает энергию системы. В дальнейшем обе площади поверхностей продолжают снижаться, хотя и менее интенсивно: ядро обращается в бислой, в котором ещё меньше места для контактов с растворителем.
Считаем меру порядка бислоя в начале и в конце моделирования.
Скачиваем индекс файл
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
%%bash
g_order -s b_md -f b_md.gpu.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
g_order -s b_md -f b_md.gpu.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')
plt.figure(figsize=(17,8))
plt.plot(np.array(start[0]), start[3], label='Start', color='red');
plt.plot(np.array(end[0]), end[3], label='End', color='green');
plt.legend();
К концу моделирования билипидный слой существенно упорядоченнее, чем вскоре после начала.